library(tidyverse)
library(lmtest)
library(sandwich)
library(texreg)
library(rstatix)
library(rsample)
library(BayesFactor)

# set working directory 
setwd("")

#### Experiment 1: IL, PA, MD 2012 ####
dCD <- read_csv("2012_congress.csv")

dCD <- mutate(dCD, vote_inc = case_when(
      House_vote_13a == 1 ~ 1,
      House_vote_13b == 1 ~ 1,
      House_vote_13c == 1 ~ 1,
      TRUE ~ 0),
      
      vote_inc_no_labels = if_else(House_vote_13a == 1, 1, 0),
      
      vote_inc_party = if_else(House_vote_13b == 1, 1, 0),
      
      vote_inc_party_incumbency = if_else(House_vote_13c == 1, 1, 0),
      
      vote_inc_no_incumbency = case_when(
            !is.na(vote_inc_no_labels) ~ vote_inc_no_labels,
            !is.na(vote_inc_party) ~ vote_inc_party,
            TRUE ~ NA_real_
      ),
      
      vote_chal_no_labels = 1 - vote_inc_no_labels,
      
      vote_chal_party = 1 - vote_inc_party,
      
      vote_chal_party_incumbency = 1 - vote_inc_party_incumbency,
      
      vote_chal_no_incumbency = 1 - vote_inc_no_incumbency,
      
      treat = case_when(
            !is.na(Chal_4A) ~ 1,
            !is.na(Chal_4B) ~ 2,
            !is.na(Attention_4C) ~ 3,
            !is.na(Attention_4D) ~ 4,
            !is.na(Attention_4E) ~ 5,
            !is.na(Attention_4F) ~ 6,
            TRUE ~ 0),
      
      treat_party = case_when(
            treat == 2 | treat == 4 ~ 1,
            treat == 0 ~ 0,
            TRUE ~ NA_real_),
      
      treat_chal = case_when(
            treat == 1 | treat == 3 ~ 1,
            treat == 0 ~ 0,
            TRUE ~ NA_real_),
      
      treat_incumb = case_when(
            treat == 5 ~ 1,
            treat == 0 ~ 0,
            TRUE ~ NA_real_),
      
      treat_chal_pt = 1 - treat_chal,
      
      chal_primed = if_else(treat != "0", 1, 0),
      
      pid2 = case_match(
            Party_ID,
            1 ~ "Democratic",
            2 ~ "Republican",
            TRUE ~ NA_character_),
      
      ind = if_else(Party_ID > 2, 1, 0),
      
      pid3 = case_match(
            Party_ID,
            1 ~ "Democratic",
            2 ~ "Republican",
            3 ~ "Independent",
            TRUE ~ NA_character_),
      
      outpartisan = case_when(
            pid2 == "Democratic" & inc_dem == 1 ~ 0,
            pid2 == "Democratic" & inc_dem == 0 ~ 1,
            pid2 == "Republican" & inc_dem == 0 ~ 0,
            pid2 == "Republican" & inc_dem == 1 ~ 1),
      
      outpartisan_subset = case_when(
            pid2 == "Democratic" & inc_dem == 1 ~ 0,
            pid2 == "Democratic" & inc_dem == 0 ~ 1,
            pid2 == "Republican" & inc_dem == 0 ~ 0,
            pid2 == "Republican" & inc_dem == 1 ~ 1,
            TRUE ~ 0),
      
      bachelors = ifelse(Education >= 5, 1, 0),
      
      state = case_when(
            supplement == 2 & cd < 7 ~ "MD",
            supplement == 2 & cd > 7 ~ "PA",
            TRUE ~ "IL"),
      
      state_cd = paste0(state, cd))

# name-only datasets
dmere <- dCD |> filter(treat %in% c(0,1,3)) |>
      mutate(treat = case_when(
            treat == 0 ~ 0,
            treat == 1 ~ 1, 
            TRUE ~ NA_real_))

dmereIL1 <- filter(dmere, state_cd == "IL1")
dmereIL3 <- filter(dmere, state_cd == "IL3")
dmereIL10 <- filter(dmere, state_cd == "IL10")
dmereMD3 <- filter(dmere, state_cd == "MD3")
dmereMD6 <- filter(dmere, state_cd == "MD6")
dmerePA10 <- filter(dmere, state_cd == "PA10")
dmerePA12 <- filter(dmere, state_cd == "PA12")

# name + party-only datasets
dnp <- dCD |> filter(treat %in% c(0,2,4)) |>
      mutate(treat = case_when(
            treat == 0 ~ 0,
            treat > 0 ~ 1, 
            TRUE ~ NA_real_))

dnpIL1 <- filter(dnp, state_cd == "IL1")
dnpIL3 <- filter(dnp, state_cd == "IL3")
dnpIL10 <- filter(dnp, state_cd == "IL10")
dnpMD3 <- filter(dnp, state_cd == "MD3")
dnpMD6 <- filter(dnp, state_cd == "MD6")
dnpPA10 <- filter(dnp, state_cd == "PA10")
dnpPA12 <- filter(dnp, state_cd == "PA12")

# name + incumbency-only datasets
dinc <- dCD |> filter(treat %in% c(0,5)) |>
      mutate(treat = case_when(
            treat == 0 ~ 0,
            treat > 0 ~ 1, 
            TRUE ~ NA_real_))

dincIL1 <- filter(dnp, state_cd == "IL1")
dincIL3 <- filter(dnp, state_cd == "IL3")
dincIL10 <- filter(dnp, state_cd == "IL10")
dincMD3 <- filter(dnp, state_cd == "MD3")
dincMD6 <- filter(dnp, state_cd == "MD6")
dincPA10 <- filter(dnp, state_cd == "PA10")
dincPA12 <- filter(dnp, state_cd == "PA12")

# for merged data 
# # name + party-only datasets
dCD4merge <- dCD |> filter(treat %in% c(0,2,4)) |>
      mutate(treat = case_when(
            treat == 0 ~ 0,
            treat > 0 ~ 1, 
            TRUE ~ NA_real_)) |> 
      select(treat, vote_chal = vote_chal_no_incumbency,
             state_cd, bachelors, outpartisan) |>
      mutate(study = state_cd) |>
      select(!state_cd)


# treat codes
# 1 = chal name
# 2 = chal name + party
# 3 = chal name + quality
# 4 = chal name + party + quality
# 5 = chal name + inc name
# 6 = chal name + inc name + party + quality

#### analysis ####

##### Name Only Treatment #####
######  Table 1 #######

## proportions hearing of challenger ##
# 3 = not heard of challenger
round(proportions(table(dmereIL1$Chal_4A)), 3)
round(proportions(table(dmereIL3$Chal_4A)), 3)
round(proportions(table(dmereIL10$Chal_4A)), 3)
round(proportions(table(dmereMD3$Chal_4A)), 3)
round(proportions(table(dmereMD6$Chal_4A)), 3)
round(proportions(table(dmerePA10$Chal_4A)), 3)
round(proportions(table(dmerePA12$Chal_4A)), 3)

# % of all respondents 
round(proportions(table(dmere$Chal_4A)), 3)

row_wise_prop_test(table(dmereIL1$vote_inc_no_labels, dmereIL1$treat_chal), detailed = T)
row_wise_prop_test(table(dmereIL3$treat_chal, dmereIL3$vote_inc_no_labels), detailed = T)
row_wise_prop_test(table(dmereIL10$treat_chal, dmereIL10$vote_inc_no_labels), detailed = T)
row_wise_prop_test(table(dmereMD3$treat_chal, dmereMD3$vote_inc_no_labels), detailed = T)
row_wise_prop_test(table(dmereMD6$treat_chal, dmereMD6$vote_inc_no_labels), detailed = T)
row_wise_prop_test(table(dmerePA10$treat_chal, dmerePA10$vote_inc_no_labels), detailed = T)
row_wise_prop_test(table(dmerePA12$treat_chal, dmerePA12$vote_inc_no_labels), detailed = T)

round(proportions(table(dmere$vote_inc_no_labels)), 3)

###### Table 2 & A1 ######
## pooled linear model ##
m1 <- lm(vote_chal_no_labels ~ treat_chal + state_cd, dCD)
texreg::screenreg(m1, digits = 3)

### bayes factor ###

dbf <- dCD |> select(vote_chal_no_labels, treat_chal, state_cd) |>
      na.omit()
mfull <- lmBF(vote_chal_no_labels ~ treat_chal + state_cd, dbf)
mnull <- lmBF(vote_chal_no_labels ~ state_cd, dbf)
mfull / mnull


##### appendix analysis #####

###### Table A3 ######
# linear model with DV that does not give party label
m <- lm(vote_chal_no_labels ~ treat_party + state_cd, dCD,
        family = binomial())
texreg::screenreg(m, digits = 3)

###### Table A4 ###### 
# linear model with only DV that gives party label
m <- lm(vote_chal_party ~ treat_party + state_cd, dCD,
        family = binomial())
texreg::screenreg(m, digits = 3)

###### Table A8 ###### 
# pooled logit #
m1 <- glm(vote_chal_no_labels ~ treat_chal + state_cd, dCD,
          family = binomial())
texreg::screenreg(m1, digits = 3)

# balance test
summary(glm(treat_chal ~ Gender + Age + factor(Race) + Education, dCD, 
            family = binomial()))

summary(glm(treat_party ~ Gender + Age + factor(Race) + Education, dCD, 
            family = binomial()))

##### Name + Party Label Treatment #####

###### Table 4 ######
## proportions hearing of challenger
round(proportions(table(dnpIL1$Chal_4B)), 3)
round(proportions(table(dnpIL3$Chal_4B)), 3)
round(proportions(table(dnpIL10$Chal_4B)), 3)
round(proportions(table(dnpMD3$Chal_4B)), 3)
round(proportions(table(dnpMD6$Chal_4B)), 3)
round(proportions(table(dnpPA10$Chal_4B)), 3)
round(proportions(table(dnpPA12$Chal_4B)), 3)


## proportions tests ##
prop_test(table(dnpIL1$treat_party, dnpIL1$vote_inc_no_incumbency), detailed = T)
prop_test(table(dnpIL3$treat_party, dnpIL3$vote_inc_no_incumbency), detailed = T)
prop_test(table(dnpIL10$treat_party, dnpIL10$vote_inc_no_incumbency), detailed = T)
prop_test(table(dnpMD3$treat_party, dnpMD3$vote_inc_no_incumbency), detailed = T)
prop_test(table(dnpMD6$treat_party, dnpMD6$vote_inc_no_incumbency), detailed = T)
prop_test(table(dnpPA10$treat_party, dnpPA10$vote_inc_no_incumbency), detailed = T)
prop_test(table(dnpPA12$treat_party, dnpPA12$vote_inc_no_incumbency), detailed = T)




#### Experiment 2: GOP Primary 2012 ####

d <- read_csv("2012_GOP.csv")

##### analysis #####

###### Table 3 ###### 
# group = 1 are rows voting for primed candidate
# n2 is group receiving prime treatment
row_wise_prop_test(table(d$vote_bachmann, d$bachmann_primed), detailed = T)
row_wise_prop_test(table(d$vote_cain, d$cain_primed), detailed = T)
row_wise_prop_test(table(d$vote_huntsman, d$huntsman_primed), detailed = T)
row_wise_prop_test(table(d$vote_perry, d$perry_primed), detailed = T)
row_wise_prop_test(table(d$vote_roemer, d$roemer_primed), detailed = T)
row_wise_prop_test(table(d$vote_santorum, d$santorum_primed), detailed = T)


### bayes factors ###
contingencyTableBF(table(d$bachmann_primed, d$vote_bachmann), "indepMulti", fixedMargin = "rows")
contingencyTableBF(table(d$vote_cain, d$cain_primed), "indepMulti", fixedMargin = "rows")
contingencyTableBF(table(d$vote_huntsman, d$huntsman_primed), "indepMulti", fixedMargin = "rows")
contingencyTableBF(table(d$vote_perry, d$perry_primed), "indepMulti", fixedMargin = "rows")
contingencyTableBF(table(d$vote_roemer, d$roemer_primed), "indepMulti", fixedMargin = "rows")
contingencyTableBF(table(d$vote_santorum, d$santorum_primed), "indepMulti", fixedMargin = "rows")


###### pooled meta-analytic estimate ######

mb <- lm(vote_bachmann ~ bachmann_primed, d)
mc <- lm(vote_cain ~ cain_primed, d)
mh <- lm(vote_huntsman ~ huntsman_primed, d)
mp <- lm(vote_perry ~ perry_primed, d)
mr <- lm(vote_roemer ~ roemer_primed, d)
ms <- lm(vote_santorum ~ santorum_primed, d)

# bootstrap 5000 samples
set.seed(02145)
samps <- rsample::bootstraps(d, 5000)

### estimate ATEs (betas) and regression standard errors (sigmas)
# bachmann betas and sigmas
lm_coefs <- function(splits, ...){
      mod <- lm(vote_bachmann ~ bachmann_primed, splits)
      coef(mod)[2]
}
lm_sigmas <- function(splits, ...){
      mod <- lm(vote_bachmann ~ bachmann_primed, splits)
      summary(mod)$sigma
}
mod_formula <- as.formula(mb)
samps$bachmann_betas <- map(.x = samps$splits, 
                            .f = lm_coefs,
                            mod_formula)
samps$bachmann_sigmas <- map(.x = samps$splits, 
                             .f = lm_sigmas,
                             mod_formula)

# cain betas and sigmas
lm_coefs <- function(splits, ...){
      mod <- lm(vote_cain ~ cain_primed, splits)
      coef(mod)[2]
}
lm_sigmas <- function(splits, ...){
      mod <- lm(vote_cain ~ cain_primed, splits)
      summary(mod)$sigma
}
mod_formula <- as.formula(mc)
samps$cain_betas <- map(.x = samps$splits, 
                        .f = lm_coefs,
                        mod_formula)
samps$cain_sigmas <- map(.x = samps$splits, 
                         .f = lm_sigmas,
                         mod_formula)

# huntsman betas and sigmas
lm_coefs <- function(splits, ...){
      mod <- lm(vote_huntsman ~ huntsman_primed, splits)
      coef(mod)[2]
}
lm_sigmas <- function(splits, ...){
      mod <- lm(vote_huntsman ~ huntsman_primed, splits)
      summary(mod)$sigma
}
mod_formula <- as.formula(mh)
samps$huntsman_betas <- map(.x = samps$splits, 
                            .f = lm_coefs,
                            mod_formula)
samps$huntsman_sigmas <- map(.x = samps$splits, 
                             .f = lm_sigmas,
                             mod_formula)


# perry betas and sigmas
lm_coefs <- function(splits, ...){
      mod <- lm(vote_perry ~ perry_primed, splits)
      coef(mod)[2]
}
lm_sigmas <- function(splits, ...){
      mod <- lm(vote_perry ~ perry_primed, splits)
      summary(mod)$sigma
}
mod_formula <- as.formula(mp)
samps$perry_betas <- map(.x = samps$splits, 
                         .f = lm_coefs,
                         mod_formula)
samps$perry_sigmas <- map(.x = samps$splits, 
                          .f = lm_sigmas,
                          mod_formula)


# roemer betas and sigmas
lm_coefs <- function(splits, ...){
      mod <- lm(vote_roemer ~ roemer_primed, splits)
      coef(mod)[2]
}
lm_sigmas <- function(splits, ...){
      mod <- lm(vote_roemer ~ roemer_primed, splits)
      summary(mod)$sigma
}
mod_formula <- as.formula(mr)
samps$roemer_betas <- map(.x = samps$splits, 
                          .f = lm_coefs,
                          mod_formula)
samps$roemer_sigmas <- map(.x = samps$splits, 
                           .f = lm_sigmas,
                           mod_formula)

# santorum betas and sigmas
lm_coefs <- function(splits, ...){
      mod <- lm(vote_santorum ~ santorum_primed, splits)
      coef(mod)[2]
}
lm_sigmas <- function(splits, ...){
      mod <- lm(vote_santorum ~ santorum_primed, splits)
      summary(mod)$sigma
}
mod_formula <- as.formula(ms)
samps$santorum_betas <- map(.x = samps$splits, 
                            .f = lm_coefs,
                            mod_formula)
samps$santorum_sigmas <- map(.x = samps$splits, 
                             .f = lm_sigmas,
                             mod_formula)

### pooled meta-analytic regression estimate ###

(mb$coefficients[2] * (1/summary(mb)$sigma) + 
            mc$coefficients[2] * (1/summary(mc)$sigma) +
            mh$coefficients[2] *  (1/summary(mh)$sigma) +
            mp$coefficients[2] * (1/summary(mp)$sigma) +
            mr$coefficients[2] * (1/summary(mr)$sigma) +
            ms$coefficients[2] * (1/summary(ms)$sigma)) / (
                  (1/summary(mb)$sigma) + 
                        (1/summary(mc)$sigma) + 
                        (1/summary(mh)$sigma) + 
                        (1/summary(mp)$sigma) + 
                        (1/summary(mr)$sigma) +
                        (1/summary(ms)$sigma)
            )

### bootstrapped meta-analytic regression standard error ###

# making data frame of betas and sigmas
samp_df <- bind_cols(
      
      bachmann_betas = list_c(samps$bachmann_betas),
      cain_betas = list_c(samps$cain_betas),
      huntsman_betas = list_c(samps$huntsman_betas),
      perry_betas = list_c(samps$perry_betas),
      roemer_betas = list_c(samps$roemer_betas),
      santorum_betas = list_c(samps$santorum_betas),
      
      bachmann_sigmas = list_c(samps$bachmann_sigmas),
      cain_sigmas = list_c(samps$cain_sigmas),
      huntsman_sigmas = list_c(samps$huntsman_sigmas),
      perry_sigmas = list_c(samps$perry_sigmas),
      roemer_sigmas = list_c(samps$roemer_sigmas),
      santorum_sigmas = list_c(samps$santorum_sigmas),
      
)

# new column is meta-analytic estimate

samp_df <- mutate(samp_df, tau = 
                        (bachmann_betas * (1/bachmann_sigmas) +
                               cain_betas * (1/cain_sigmas) +
                               huntsman_betas * (1/huntsman_sigmas) +
                               perry_betas * (1/perry_sigmas) +
                               roemer_betas * (1/roemer_sigmas) +
                               santorum_betas * (1/santorum_sigmas)) / (
                                     (1/bachmann_sigmas) + 
                                           (1/cain_sigmas) + 
                                           (1/huntsman_sigmas) + 
                                           (1/perry_sigmas) +
                                           (1/roemer_sigmas) + 
                                           (1/santorum_sigmas)
                               ))

# bootstrapped standard error     
sd(samp_df$tau)

# bootstrapped 95% CI
quantile(samp_df$tau, .025)
quantile(samp_df$tau, .975)

##### appendix analysis #####

###### controling for hispanic ######
mb <- lm(vote_bachmann ~ bachmann_primed + hispanic, d)
mc <- lm(vote_cain ~ cain_primed + hispanic, d)
mh <- lm(vote_huntsman ~ huntsman_primed + hispanic, d)
mp <- lm(vote_perry ~ perry_primed + hispanic, d)
mr <- lm(vote_roemer ~ roemer_primed + hispanic, d)
ms <- lm(vote_santorum ~ santorum_primed + hispanic, d)


# bootstrap 5000 samples
samps <- rsample::bootstraps(d, 5000)

### estimate ATEs (betas) and regression standard errors (sigmas)
# bachmann betas and sigmas
lm_coefs <- function(splits, ...){
      mod <- lm(vote_bachmann ~ bachmann_primed + hispanic, splits)
      coef(mod)[2]
}
lm_sigmas <- function(splits, ...){
      mod <- lm(vote_bachmann ~ bachmann_primed + hispanic, splits)
      summary(mod)$sigma
}
mod_formula <- as.formula(mb)
samps$bachmann_betas <- map(.x = samps$splits, 
                            .f = lm_coefs,
                            mod_formula)
samps$bachmann_sigmas <- map(.x = samps$splits, 
                             .f = lm_sigmas,
                             mod_formula)

# cain betas and sigmas
lm_coefs <- function(splits, ...){
      mod <- lm(vote_cain ~ cain_primed + hispanic, splits)
      coef(mod)[2]
}
lm_sigmas <- function(splits, ...){
      mod <- lm(vote_cain ~ cain_primed + hispanic, splits)
      summary(mod)$sigma
}
mod_formula <- as.formula(mc)
samps$cain_betas <- map(.x = samps$splits, 
                        .f = lm_coefs,
                        mod_formula)
samps$cain_sigmas <- map(.x = samps$splits, 
                         .f = lm_sigmas,
                         mod_formula)

# huntsman betas and sigmas
lm_coefs <- function(splits, ...){
      mod <- lm(vote_huntsman ~ huntsman_primed + hispanic, splits)
      coef(mod)[2]
}
lm_sigmas <- function(splits, ...){
      mod <- lm(vote_huntsman ~ huntsman_primed + hispanic, splits)
      summary(mod)$sigma
}
mod_formula <- as.formula(mh)
samps$huntsman_betas <- map(.x = samps$splits, 
                            .f = lm_coefs,
                            mod_formula)
samps$huntsman_sigmas <- map(.x = samps$splits, 
                             .f = lm_sigmas,
                             mod_formula)


# perry betas and sigmas
lm_coefs <- function(splits, ...){
      mod <- lm(vote_perry ~ perry_primed + hispanic, splits)
      coef(mod)[2]
}
lm_sigmas <- function(splits, ...){
      mod <- lm(vote_perry ~ perry_primed + hispanic, splits)
      summary(mod)$sigma
}
mod_formula <- as.formula(mp)
samps$perry_betas <- map(.x = samps$splits, 
                         .f = lm_coefs,
                         mod_formula)
samps$perry_sigmas <- map(.x = samps$splits, 
                          .f = lm_sigmas,
                          mod_formula)


# roemer betas and sigmas
lm_coefs <- function(splits, ...){
      mod <- lm(vote_roemer ~ roemer_primed + hispanic, splits)
      coef(mod)[2]
}
lm_sigmas <- function(splits, ...){
      mod <- lm(vote_roemer ~ roemer_primed + hispanic, splits)
      summary(mod)$sigma
}
mod_formula <- as.formula(mr)
samps$roemer_betas <- map(.x = samps$splits, 
                          .f = lm_coefs,
                          mod_formula)
samps$roemer_sigmas <- map(.x = samps$splits, 
                           .f = lm_sigmas,
                           mod_formula)

# santorum betas and sigmas
lm_coefs <- function(splits, ...){
      mod <- lm(vote_santorum ~ santorum_primed + hispanic, splits)
      coef(mod)[2]
}
lm_sigmas <- function(splits, ...){
      mod <- lm(vote_santorum ~ santorum_primed + hispanic, splits)
      summary(mod)$sigma
}
mod_formula <- as.formula(ms)
samps$santorum_betas <- map(.x = samps$splits, 
                            .f = lm_coefs,
                            mod_formula)
samps$santorum_sigmas <- map(.x = samps$splits, 
                             .f = lm_sigmas,
                             mod_formula)

### pooled meta-analytic regression estimate ###

(mb$coefficients[2] * (1/summary(mb)$sigma) + 
            mc$coefficients[2] * (1/summary(mc)$sigma) +
            mh$coefficients[2] *  (1/summary(mh)$sigma) +
            mp$coefficients[2] * (1/summary(mp)$sigma) +
            mr$coefficients[2] * (1/summary(mr)$sigma) +
            ms$coefficients[2] * (1/summary(ms)$sigma)) / (
                  (1/summary(mb)$sigma) + 
                        (1/summary(mc)$sigma) + 
                        (1/summary(mh)$sigma) + 
                        (1/summary(mp)$sigma) + 
                        (1/summary(mr)$sigma) +
                        (1/summary(ms)$sigma)
            )

### bootstrapped meta-analytic regression standard error ###

# making data frame of betas and sigmas
samp_df <- bind_cols(
      
      bachmann_betas = list_c(samps$bachmann_betas),
      cain_betas = list_c(samps$cain_betas),
      huntsman_betas = list_c(samps$huntsman_betas),
      perry_betas = list_c(samps$perry_betas),
      roemer_betas = list_c(samps$roemer_betas),
      santorum_betas = list_c(samps$santorum_betas),
      
      bachmann_sigmas = list_c(samps$bachmann_sigmas),
      cain_sigmas = list_c(samps$cain_sigmas),
      huntsman_sigmas = list_c(samps$huntsman_sigmas),
      perry_sigmas = list_c(samps$perry_sigmas),
      roemer_sigmas = list_c(samps$roemer_sigmas),
      santorum_sigmas = list_c(samps$santorum_sigmas),
      
)

# new column is meta-analytic estimate

samp_df <- mutate(samp_df, tau = 
                        (bachmann_betas * (1/bachmann_sigmas) +
                               cain_betas * (1/cain_sigmas) +
                               huntsman_betas * (1/huntsman_sigmas) +
                               perry_betas * (1/perry_sigmas) +
                               roemer_betas * (1/roemer_sigmas) +
                               santorum_betas * (1/santorum_sigmas)) / (
                                     (1/bachmann_sigmas) + 
                                           (1/cain_sigmas) + 
                                           (1/huntsman_sigmas) + 
                                           (1/perry_sigmas) +
                                           (1/roemer_sigmas) + 
                                           (1/santorum_sigmas)
                               ))

# bootstrapped standard error     
sd(samp_df$tau)

# bootstrapped 95% CI
quantile(samp_df$tau, .025)
quantile(samp_df$tau, .975)



#### Experiment 3: MA 2018 ####

#setwd("Dropbox/Name recognition/")

dMA <- read_csv("MA_SOS_2018.csv")

dMA <- mutate(dMA, vote_dem = case_match(
      q5,
      1 ~ 1, 
      2 ~ 0, 
      3 ~ NA_real_),
      
      vote_chal = 1 - vote_dem,
      
      treat2 = case_match(
            version,
            1 ~ "Dem",
            2 ~ "GOP",
            3 ~ "Dem + Party",
            4 ~ "GOP + Party",
            5 ~ "Control"),
      
      treat = case_match(version,
                         4 ~ 1,
                         5 ~ 0,
                         TRUE ~ NA_real_),
      
      pid3 = case_match(
            PID,
            1 ~ "Democratic",
            2 ~ "Republican",
            3 ~ "Independent",
            TRUE ~ NA_character_),
      
      bachelors = case_match(q10,
                             c(5,6) ~ 1, 
                             .default = 0),
      
      outpartisan = case_match(
            PID,
            1 ~ 0,
            2 ~ 1,
            TRUE ~ NA_real_),
      
      pid7ish = case_when(
            PID == 1 & q3 == 2 ~ 1,
            PID == 1 & q3 == 1 ~ 2,
            PID == 3 & q3 == 2 ~ 3, # ind + disapprove of Trump
            PID == 3 & q3 == 1 ~ 4, # ind + approve of Trump
            PID == 2 & q3 == 2 ~ 5,
            PID == 2 & q3 == 1 ~ 6,
            TRUE ~ NA_real_)
)

#### analysis ####

##### Name + Party Treatment #####

###### Table 4 cont. ###### 
## proportions tests ##
row_wise_prop_test(table(dMA$vote_chal, dMA$treat), 
                   detailed = T, p.adjust.method = "none")


### for merged dataset
dMA4merge <- dMA |> filter(version %in% c(4, 5)) |> 
      select(treat, vote_chal, outpartisan, bachelors) |>
      mutate(study = "MA 2018")


#### Experiment 4: KY 2011 ####
dKY <- read_csv("2011_KY.csv")

dKY <- mutate(dKY, vote_ppool = case_match(
      vote_intent,
      1 ~ 0, 
      2 ~ 0,
      3 ~ 1, 
      4 ~ 1, 
      5 ~ NA_real_), # 1 and 2 are voting for incumbent
      
      vote_inc = ifelse(vote_intent < 3, 1 ,0),
      
      vote_chal = 1 - vote_inc,
      
      pid7 = na_if(PID, 8),
      
      pid3 = case_when(
            pid7 < 4 ~ "Democratic",
            pid7 == 4 ~ "Independent",
            pid7 > 4 ~ "Republican",
            pid7 == 8 ~ NA_character_), 
      
      bachelors = case_match(educ,
                             c(6,7) ~ 1,
                             8 ~ NA_real_, 
                             .default = 0),
      
      treat = ifelse(treatment_favorability == 1, 1, 0),
      
      outpartisan = case_when(
            pid7 < 4 ~ 0,
            pid7 == 4 ~ NA_real_,
            pid7 > 4 ~ 1,
            pid7 == 8 ~ NA_real_))

dKY4merge <- dKY |> select(treat, vote_chal, bachelors, outpartisan) |>
      mutate(study = "KY 2011")

##### analysis #####

###### Table 4 cont. ###### 
## proportions tests ##
# 64% never heard of p'pool
prop.table(table(dKY$ppool_favorability))

proportions(table(dKY$vote_chal, dKY$treat), 2)

row_wise_prop_test(table(dKY$vote_chal, dKY$treat), 
                   p.adjust.method = "none", detailed = T)


#### merge datasets ####

dbig <- bind_rows(dCD4merge, dMA4merge, dKY4merge)

##### pooled proportions test #####

row_wise_prop_test(table(dbig$vote_chal, dbig$treat), 
                   p.adjust.method = "none", detailed = T)

###### Table 5 & A2 ######
m5 <- lm(vote_chal ~ treat + study, dbig)
texreg::screenreg(m5, digits = 3)

### bayes factor ###

dbf <- dbig |> select(vote_chal, treat, study) |>
      na.omit()
mfull <- lmBF(vote_chal ~ treat + study, dbf)
mnull <- lmBF(vote_chal ~ study, dbf)
mfull / mnull


###### Table A5 ###### 
# interaction with Out-Party ID 
m5 <- glm(vote_chal ~ treat*outpartisan + study, dbig,
          family = binomial())
texreg::screenreg(m5, digits = 3)

###### Table A6 ###### 
# interaction with education
m5 <- glm(vote_chal ~ treat*bachelors + study, dbig,
          family = binomial())
texreg::screenreg(m5, digits = 3)

###### Table A7 ###### 
m5 <- glm(vote_chal ~ treat*study, dbig,
          family = binomial())
texreg::screenreg(m5, digits = 3)

###### Table A9 ######

m5 <- glm(vote_chal ~ treat + study, dbig,
          family = binomial())
texreg::screenreg(m5, digits = 3)

###### power analysis ######

# STATA command for merged dataset and halved effect size
# power twoproportions .36 .4250, test(chi2) n1(2413) n2(2204)



